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Abstract 

Background: The molecular mechanisms involved in plant tolerance to either drought or cold have been extensively 
studied in many plant species. However, few studies have focused on their comparisons especially using non-model 
plants with strong tolerance to both stresses. Ammopiptanthus mongolicus (Maxim, ex Kom.) Cheng f. is the only 
evergreen broadleaf shrub grown in the central Asian desert and it has very strong cold and drought tolerance. To 
provide further insights into plant tolerance, the transcriptome profiles of drought- and cold-treated A mongolicus 
seedlings were analyzed using lllumina technology and differentially expressed genes (DEGs) were compared. 

Results: A comprehensive transcriptome of A. mongolicus was sequenced using pooled mRNA extracted from 
drought-, cold-stressed and unstressed seedlings as well as leaves from naturally grown shrub. These sequences were 
assembled into 86058 unigenes, of which 51014 unigenes had an annotated function and 2440 encoded transcription 
factors (TPs). Transcriptome profiles were analyzed in A. mongolicus seedlings after drought and cold treatments at 
three time points (2, 8 and 24 h). Between 3917 and 6102 unigenes were identified as DEGs at a single time point in 
both stresses. Among these DEGs 2028 and 2026 DEGs were common across the three time points of drought and 
cold treatments respectively, and 971 DEGs were co-regulated by both stresses. Functional enrichment analyses 
identified many common or specific biological processes and gene sets in response to drought and cold stresses. The 
most pronounced findings are that flavonoid biosynthesis genes were enriched in the DEGs co-up-regulated by both 
stresses; while membrane protein genes and genes related to chloroplast were abundant in the DEGs specifically 
up-regulated by drought or cold, respectively. Furthermore, the DREB, ERF, NAG and WRKYTFs were predominantly 
co-up-regulated by both stresses. 

Conclusions: The present study provides the most comprehensive transcriptome resource and the first dynamic 
transcriptome profiles of A. mongolicus under drought and cold stresses. This information will deepen our 
understanding of plant tolerance to drought and cold. The up-regulated DEGs will be valuable for further investigations 
of key genes and molecular mechanisms involved in the adaptation of A mongolicus to harsh environments. 

Keywords: Ammopiptanthus mongolicus, Drought, Cold, Transcriptome profiling, lllumina sequencing 



^BMC 

Genomics 



* Correspondence: xinguo.li(S)csiro.au; wangmaoyan(S) 163.com 
^Equal contributors 

^Agriculture Flagship, The Commonwealth Scientific and Industrial Research 
Organization (CSIRO), Canberra ACT 2601, Australia 
^College of Life Sciences, Inner Mongolia Agricultural University, No. 306 
Zhaowuda Street, Hohhot 010018, China 

O© 2014 Wu et a!.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative 
BIoIVIGCI Central commons Attribution License (http://creativecommons.Org/licenses/by/2.0), which permits unrestricted use, distribution, and 
reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain 
Dedication waiver (http://creativecommons.0rg/publicdomain/zero/l.O/) applies to the data made available in this article, 
unless otherwise stated. 



Wu et al. BMC Genomics 2014, 15:671 
http://www.bionnedcentral.conn/1471 -21 64/1 5/671 



Page 2 of 16 



Background 

Drought and low temperature are two major environmen- 
tal stresses that greatly affect plant growth and crop pro- 
duction worldwide. This is particularly the case in arid 
and cold areas. To improve the tolerance of crops and 
other economic plants to these stresses, it is essential to 
know how plants respond to these stresses and which 
genes and biological pathways are involved in the stress 
tolerance. Over the past decades, thousands of genes and 
dozens of metabolic and signaling pathways have been 
identified in plants during drought and cold stresses [1-4]. 
Some of these genes have been confirmed to significantly 
improve plant tolerance to drought, cold and/or salt 
stresses [5,6]. In addition, many genes have been found to 
be commonly or specifically regulated under drought and 
cold in Arabidopsis [7], rice [8] and Brachypodium [9] at 
transcriptional level. However, similar studies are very lim- 
ited in other plants, and the molecular mechanisms con- 
trolling plant tolerance to both drought and cold remain 
largely unknown. 

Ammopiptanthus is an endangered survivor from the 
Tethys in the Tertiary Period, and is the only evergreen 
broadleaf shrub grown in the central Asian desert. The 
genus Ammopiptanthus (Leguminosae) comprises of two 
species: A, mongolicus (Maxim, ex Kom.) Cheng f. and A, 
nanus (M. Pop.) Cheng f. Both species are diploid with 18 
chromosomes and have high stress tolerance [10]. The 
habitats of A, mongolicus are marked by arid climate with 
an annual precipitation less than 200 mm but a mean an- 
nual evaporation up to 3000 mm. Its natural distribution 
areas are also characterized as sandy or stony soil with 
high salinity, intense ultraviolet irradiation, and seasonally 
extreme temperature from about -30°C in winter to more 
than 40°C in summer. The extreme tolerance of this spe- 
cies to harsh environments makes it invaluable for explor- 
ing key stress-tolerant genes and mechanisms, especially 
those involved in cold and drought tolerance. 

The transcriptomes of A, mongolicus roots (treated by 
20% PEG-6000 for 72 hours) and seedlings (cultured at 
normal condition or 4°C for 14 days) were recently se- 
quenced by 454 pyrosequencing [11] and Illumina tech- 
nology [12], respectively. These studies generated 29056 
and 82795 unigenes and identified 32728 cold-regulated 
genes. However, transcriptome profiling of this species 
under drought stress and its comparison with cold stress 
have not been reported in the literatures. 

In the present study, we firstly performed a comprehen- 
sive transcriptome sequencing of A, mongolicus using 
pooled mRNA extracted from drought-treated, cold-treated 
and control seedlings as well as young leaves of plants 
grown in desert in both summer and winter. We then ap- 
plied RNA-Seq (RNA sequencing) to investigate seven 
cDNA libraries derived from the seedling samples exposed 
to drought and cold stresses for 2 h, 8 h and 24 h, and the 



non-stressed control seedlings respectively. Finally, we ana- 
lyzed the DEGs and identified common and specific DEGs 
during drought and cold treatments. Our study provided 
the first transcriptome profile of A, mongolicus under 
drought stress and the comparison of it with the transcrip- 
tome profile under cold stress. 

Results 

Sequencing, de novo assembly and functional annotation 
of the A, mongolicus transcriptome 

Illumina sequencing of the pooled mRNA generated 47.29 
million (M) clean paired-end reads with Q20 over 98%. De 
novo assembly of the clean reads resulted in 86058 uni- 
genes with an average length of 756 bp and N50 length of 
1279 bp (Table 1). Approximately 41.4% and 22.6% of 
these unigenes had a length > 500 bp or > 1000 bp, re- 
spectively (Figure 1). 

Of the 86058 sequences, 51014 (59.3%) unigenes had at 
least one significant match in the NCBI protein databases 
(E-value < le-5). Distribution of the annotated unigenes in 
each database was shown in the Additional file 1. By search- 
ing against the PlantTFDB (http://planttfdb.cbi.edu.cn) 
(E-value < le-5) in combination with annotations in the 
NR (non-redundant) and Swiss-Prot databases, 2440 
unigenes were annotated as transcription factors (TPs), 
representing 4.8% of the annotated unigenes and covering 
57 TP families (Table 2). Furthermore, 47486 unigenes 
had significant hits in the TIGR Plant Transcript Assem- 
blies database (E-value < le-5). Most of these hits were 
from Glycine max, one of its closest relatives, followed by 
Medicago truncatula, Lotus japonicus, Malus x, Glycine 
soja, and others. Surprisingly, only 919 (1.9%) and 787 
(1.7%) unigenes had significant similarity to the sequences 
of Arabidopsis thaliana and Oryza sativa, respectively 
(Additional file 2). 

Of the 86058 unigenes, 23834 unigenes were assigned 
to 44 GO (Gene Ontology) terms. The terms "Cell", 
"Cell part", "Organelle", "Binding", "Catalytic activity", 
"Metabolic process", "Cellular process" and "Response to 
stimulus" were mostly dominant (Figure 2). Moreover, 

Table 1 Overview of the transcriptome sequencing 
and assembly 



Item 

Total number of raw reads 
Total number of clean reads 
Total clean nucleotides (bp) 
Q20 percentage (%) 
Percentage (%) of assembled reads 
Number of contigs or unigenes 
Average length (bp) 
N50 (bp) 



Reads Contigs Unigenes 

50477558 

47291846 

4256266140 - 

98.36 

99.84 

190185 86058 
325 756 
525 1279 
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Figure 1 Size distribution of the unigene sequences. The 

numbers of six groups of unigenes with different lengths (200 to 
499 bp, 500 to 999 bp, 1000 to 1499 bp, 1500 to 1999 bp, 2000 to 
2999 bp and > 3000 bp) were shown, respectively. The percentages 
of the unigenes in each group out of the total unigenes (86058) 
were also indicated. 



17100 unigenes were classified into 25 COG (Clusters of 
Orthologous Groups) classes. The class "General function 
prediction only" represented the largest group, followed 
by the classes "Transcription", "Signal transduction mech- 
anisms" and "Amino acid transport and metabolism", etc. 
Notably, nearly 1600 of the annotated unigenes were 
assigned to the class of "Function unknown" (Additional 
file 3). 

To gain insights into function interactions among the 
unigenes, 26999 of the 86058 unigenes were assigned 
to 125 KEGG (Kyoto Encyclopedia of Genes and 
Genomes) pathways. Of the top 10 pathways (Table 3), 
"Metabolic pathways" were mostly represented, followed 
by "Biosynthesis of secondary metabolites", "Plant- 
pathogen interaction" and "Plant hormone signal trans- 
duction", etc, Unigenes encoding core components in 
the plant hormone signal transduction pathways were 
present, including ABA (abscisic acid), BR (brassinoster- 
oid) and ethylene signaling pathways. These pathways 
were well-documented to participate in the response 
and adaptation of plants to stressful environments 
[13,14]. 

Identification and validation of differentially expressed 
genes (DEGs) 

Transcriptome changes of the A. mongolicus seedlings 
under drought and cold stresses (2 h, 8 h and 24 h) and 



Table 2 Transcription factor (TF) gene families identified in the A, mongolicus transcriptome^ 



TF family 


Number of unigenes 


TF family 


Number of unigenes 


TF family 


Number of unigenes 


bHLH 


243 


GTE 


39 


Nin-like 


9 


MYB 


196 


HSE 


38 


HB-PHD 


8 


AP2/EREBP 


172 


Dof 


35 


WOX 


7 


C2H2 


156 


JUMONJI 


33 


YABBY 


7 


MYB-related 


155 


TCP 


32 


BESl 


6 


WRKY 


149 


LBD 


27 


Whirly 


6 


bZIP 


127 


NE-YC 


25 


GeBP 


4 


C3H 


82 


SBP 


24 


SRS 


4 


GRAS 


74 


NE-YA 


18 


PLATZ 


4 


G2-like 


68 


ZE-HD 


18 


BBR/BPC 


3 


NAC 


60 


GRE 


17 


EIL 


3 


CAMTA 


59 


co-like 


16 


SlEa-like 


3 


PARI 


57 


NE-YB 


16 


VOZ 


3 


TALE 


57 


DBB 


15 


LIM 


2 


HD-ZIP 


56 


E2E-DP 


15 


NE-Xl 


1 


ARE 


54 


M-type 


14 


STAT 


1 


Trihelix 


52 


MIKC 


12 


KAN 


1 


B3/ABI3 


43 


ARR-B 


11 






GATA 


42 


TUB 


11 






HB-other 


40 


CPP 


10 






*A total of 2440 TF unigenes were identified, covering 57 TF families. 
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Figure 2 Histogram of GO classification of tlie annotated unigenes. Unigenes with best BLAST hits were aligned to GO database. In total, 
23834 unigenes were assigned to at least one GO term and were grouped into three main GO categories and 44 GO terms. Left Y-axis 
represented the percentages of unigenes in each main category. Right Y-axis indicated the numbers of unigenes in each GO term. 



unstressed control were analyzed using Illumina technol- 
ogy. Between 8.29 and 8.87 M clean reads (50 bp in length) 
were obtained from the seven libraries. Approximately 7.18 
to 7.78 M of the clean reads in each library were mapped 
to the 86058 unigenes assembled earlier in this study (mis- 
match < 2 bp), accounting for 85.0% to 88.4% of the total 
clean reads (Table 4). In the Dl, D2 and D3 libraries (2, 8 
and 24 h after drought treatment), 2565, 2703 and 2589 
unigenes were up-regulated and 2193, 3160 and 3513 uni- 
genes were down-regulated, respectively. Similar results 
were observed in the CI, C2 and C3 libraries (2, 8 and 24 h 
after cold treatment), including 2027, 1894 and 1924 uni- 
genes with up-regulation and 2011, 2023 and 2489 uni- 
genes with down-regulation, respectively (Figure 3). Taken 



together, 9487 and 6524 unigenes were identified as DEGs 
with an absolute value of log2 Ratio > 1 and FDR (false dis- 
covery rate) < 0.001 at one or more time points during 
drought and cold stresses, respectively. 

To validate the RNA-Seq results, 11 DEGs with different 
expression patterns were selected for RT-qPCR (real-time 
quantitative PGR) analysis. These genes encode NAG 
(NAM, ATAF and GUG) and DREB (dehydration respon- 
sive element binding) TFs, HSP (heat shock protein), GLP 
(germin-like protein), TPP (trehalose-6-phosphate phos- 
phatase) and GGD (carotenoid cleavage dioxygenase). 
Nine (81.8%) of these genes showed similar expression 
patterns to those detected by the RNA-Seq data (Figure 4), 
while other two genes had similar expression trends at 



Table 3 The top 10 enriched pathways in the A. mongolicus transcriptome^ 



Pathway 


Number of unigenes 


Percentage (%) 


Pathway ID 


Metabolic pathways 


5623 


20.8 


koOllOO 


Biosynthesis of secondary metabolites 


2723 


10.1 


kolllO 


Plant-pathogen interaction 


1753 


6.5 


ko04626 


Plant hormone signal transduction 


1740 


6.4 


ko04075 


Purine metabolism 


1483 


5.5 


ko00230 


RNA transport 


1124 


4.2 


ko03013 


Endocytosis 


1004 


3.7 


ko04144 


Splicesome 


945 


3.5 


ko03040 


Ribosome 


908 


3.4 


ko03010 


Glycerophospholipid metabolism 


904 


3.4 


ko00564 



^Percentage (%): number of the unigenes assigned to certain pathway/total 26999 unigenes assigned to 125 pathways. 
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Table 4 Statistics of seven libraries and their reads mapping^ 



Library 


CK 


D1 


D2 


D3 


CI 


C2 


C3 


Raw reads (M) 


8.91 


8.34 


8.66 


8.92 


8.59 


8.38 


8.37 


Clean reads (M) 


8.86 


8.29 


8.62 


8.87 


8.52 


8.32 


8.32 


Mapped reads (M) 


7.53 


7.31 


7.62 


7.78 


7.42 


7.18 


7.23 


% of mapped reads 


85.0 


88.2 


88.4 


87.7 


87.1 


86.3 


86.9 


Perfectly matched reads (M) 


6.25 


6.12 


6.37 


6.53 


6.22 


5.90 


5.98 


% of perfectly matched reads 


71.6 


73.8 


73.9 


74.2 


72.9 


70.6 


71.9 


Mismatched reads (M) 


1.28 


1.19 


1.25 


1.25 


1.20 


1.29 


1.24 


% of mismatched reads 


14.7 


14.4 


14.5 


14.2 


14.1 


15.5 


15.0 


Uniquely matched reads (M) 


5.60 


5.27 


5.54 


5.75 


5.46 


5.27 


5.31 


% of uniquely matched reads 


63.2 


63.6 


64.3 


64.8 


64.1 


63.3 


63.8 


Unmapped reads (M) 


1.2 


0.98 


1.00 


1.02 


1.10 


1.14 


1.09 


% of unmapped reads 


13.5 


11.8 


11.6 


11.5 


12.9 


13.7 


13.1 


Total unigenes 


60279 


61029 


60149 


59789 


61231 


62693 


61568 


Coverage (%) 


70.0 


70.9 


69.9 


69.5 


71.2 


72.9 


71.5 



*CK: untreated control; D1, D2 and D3: 2, 8 and 24 h after drought treatment; CI, C2 and C3: 2, 8 and 24 h after cold treatment; M: million; Mismatch: < 2 bp; 
Coverage (%): total unigenes in each library/86058 assembled unigenes. 



one to two time points. Thus, the RNA-Seq results were 
considerably reliable for the identification of DEGs during 
drought and cold stresses in this study. 

Differentially expressed genes in response to drought 
stress 

Of 4758 to 6102 DEGs identified at the three time points of 
drought treatment, 2028 DEGs (21.4% of total DEGs) were 
common across all three time points, including 779 DEGs 
with up-regulation, 1185 DEGs with down-regulation 
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Figure 3 The numbers of DEGs up- and down-regulated in each 
library. The DEGs (differentially expressed genes) were identified using 
a threshold of FDR (false discovery rate) < 0.001 and absolute value of 
log2 Ratio > 1. Dl, D2 and D3 represented three drought-treated 
libraries, corresponding to 2, 8 and 24 h after drought treatment; CI, 
C2 and C3 represented three cold-treated libraries, corresponding to 2, 
8 and 24 h after cold treatment. The numbers of the DEGs in each 
library were shown in terms of up-regulation and down-regulation, 
respectively. 



(Figure 5) and 64 DEGs with complex regulation. The com- 
mon DEGs represented reproducible changes at each time 
point Between 1218 and 1874 DEGs (12.8% to 19.8% of 
total DEGs) were specifically regulated at a single time 
point, including 600 to 1178 DEGs with up-regulation and 
591 to 1038 DEGs with down-regulation. Besides, many 
DEGs were common at two of the three time points 
(Figure 5). The top 100 DEGs commonly regulated by 
drought stress were listed in the Additional file 4. 

Among the common DEGs up-regulated at the three 
time points, genes known to respond to drought or other 
stresses were largely represented, such as those encoding 
ripening-related proteins, LEA (late embryogenesis abun- 
dant) proteins, peroxidases, transporters, enzymes in 
the flavonoid biosynthetic pathways, protein kinases, 
and ethylene receptors. Many TFs in multiple fam- 
ilies were also up-regulated, including the AP2/EREBP 
(APETALA2/ethylene-responsive element binding pro- 
tein), NAG, WRKY, MYB (myeloblastosis oncogene), 
bHLH (basic helix-loop-helix), C2H2 (Cys2/His2-type 
zinc finger protein 2), bZIP (basic-domain leucine- 
zipper protein), C3H (Cys3/His-type zinc finger pro- 
tein), GRAS (GAT, RGA and SCR), TCP (TBI, CYC and 
PCF), COL (constans-like), Dof (DNA binding with one 
finger), and GRF (general regulatory factor) families, etc. 
Approximately 17.2% (134) of the common DEGs were 
functional unknowns. Similarly, the common DEGs 
down-regulated at the three time points of drought 
stress had a wide range of functions. For example, genes 
encoding aquaporin NIP, polyol transporter, gibberellin, 
small heat shock proteins, calcium-binding proteins, 
and some bHLH TFs appeared in this group. Interest- 
ingly, the unknowns in the down-regulated DEGs (397, 
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Gene name 

Figure 4 Validation of expression patterns of nine DEGs by RT-qPCR assay. Nine DEGs (differentially expressed genes) showed similar 
expression patterns between RNA-Seq data (A) and RT-qPCR assay (B). These genes included six genes up-regulated under drought and cold 
stresses: AmNACl (NAM, ATAF and CUC), AmNAC2, AmNACl 1, AmDREB2C (dehydration responsive element binding 2C), AmDREB2 and AmTPP 
(trehalose-6-phosphate phosphatase, not significantly down-regulated in CI library); two genes down-regulated under both stresses: AmtiSP (heat 
shock protein) and AmCCD (carotenoid cleavage dioxygenase); and one gene with different expression trends: AmGLP (germin-like protein). Dl, 
D2 and D3 represented 2, 8 and 24 h after drought treatment; CI, C2 and C3 represented 2, 8 and 24 h after cold treatment. The total RNAs used 
for the transcriptome profile analyses using RNA-Seq were used in the RT-qPCR assay, including three biological replicates. The relative expression 
levels of the selected genes were calculated using the 2"^^^ method. Error bars represented the standard deviation of the mean expression 
values. 



33.5%) were almost twice as those in the commonly up- 
regulated DEGs. 

Based on GO enrichment analysis, a large number of the 
DEGs commonly up-regulated in drought were mostly 
enriched in two similar biological processes: "Response to 
stimulus" and "Response to stress", followed by "Organic 
substance transport", "Cellular aromatic compound meta- 
bolic process" and "Protein phosphorylation". In contrast, 
the commonly down-regulated DEGs were significantly 
enriched in the biological processes "Monocarboxylic acid 
metabolism", "Small molecule biosynthesis" and 



"Nucleobase transport". In the cellular component and mo- 
lecular function categories, the commonly up- and down- 
regulated DEGs were significantly enriched in similar 
terms, such as "Membrane" and "Oxidoreductase activity" 
(Additional file 5). KEGG pathway enrichment analysis re- 
vealed that genes involved in microbial metabolism in di- 
verse environments and some secondary metabolite 
degradation were significantly enriched in both the up- and 
down-regulated DEGs. However, genes in the pathways for 
nitrogen and phenylalanine metabolisms were specifically 
enriched in the up-regulated DEGs; while genes involved in 
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Figure 5 Venn diagrams showing the numbers of common and specific DEGs at different time points of drought and cold treatments. The 

up-regulated (A) and down-regulated (B) DEGs (differentially expressed genes) identified at 2, 8 and 24 h after drought (Dl, D2 and D3) and cold 
(CI, C2 and C3) treatments were analyzed. The numbers of common and specific DEGs at different time points were shown in the overlapping and 
non-overlapping regions, respectively. The total numbers of up- and down-regulated DEGs at each time point were indicated in parentheses. 

V J 



metabolic pathways and the pathways for starch and sucrose 
metabolism as well as carotenoid biosynthesis were only 
enriched in the down-regulated DEGs (Additional file 6). 

As expected, many regulatory genes were specifically 
up-regulated at the first time point (2 h), such as genes en- 
coding some types of protein kinases, calcium-binding 
proteins, TPs, as well as enzymes or components in plant 
hormone (especially ABA) biosyntheses or signaling path- 
ways. After the treatment was prolonged (8 h to 24 h), 
many of these regulatory genes were no longer differen- 
tially expressed. Instead, a large number of functional 
genes were up-regulated, such as genes encoding osmo- 
protectant synthetases, dehydration-protective proteins, 
PR (pathogenesis-related) proteins, hydrolases, and de- 
toxification enzymes. 

Differentially expressed genes in response to cold stress 

Between 3917 and 4413 DEGs were identified at a single 
time point during cold treatment, and 2026 DEGs were 
common across the three time points, including 851 with 
up-regulation and 1174 with down-regulation (Figure 5). 
The top 100 DEGs commonly regulated by cold stress were 
listed in the Additional file 7. Similar to drought stress, 
many commonly up-regulated DEGs were annotated to 
well-known stress- or cold-induced genes. Some of these 
genes are identical to the drought-induced genes; while 
others encode different types of proteins, such as fatty acid 



desaturases, chitinases and proteins or enzymes located in 
chloroplast. The common DEGs down-regulated at three 
time points of cold stress also encode various fianctional 
and regulatory proteins. Notably, 16.8% (143) and 40.3% 
(473) of the commonly up- or down-regulated DEGs were 
unknowns, respectively. 

A large number of the DEGs up-regulated under cold 
stress were significantly enriched in the biological pro- 
cesses "Response to stimulus", "Response to stress", 
"Defense response" and "Cellular aromatic compound 
metabolic process". In the down-regulated DEGs, the 
mostly represented GO term of biological process was 
also "Response to stimulus", followed by "Lipid biosyn- 
thetic process", "Organic acid biosynthetic process" and 
"Lipid localization". Interestingly, in the cellular compo- 
nent category, the up-regulated DEGs were mostly 
enriched in "Plastid" and "Chloroplast"; while the down- 
regulated DEGs were not significantly enriched in any 
terms of this category (Additional file 8). KEGG path- 
way analysis revealed that the metabolic pathways and 
the pathways for secondary metabolite biosynthesis or 
degradation were significantly enriched in both the up- 
and down-regulated DEGs. However, the pathways in- 
volved in phenylalanine metabolism and photosynthesis 
were exclusively enriched in the up-regulated DEGs. In 
contrast, the pathways related to fatty acid metabolism, 
carotenoid biosynthesis and several amino acid 
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degradations were only enriched in the down-regulated 
DEGs (Additional file 9). 

Differentially expressed genes co-regulated during 
drought and cold stresses 

Among the 9487 drought-regulated and 6524 cold- 
regulated DEGs with significant expression changes 
(absolute value of log2 Ratio > 1 and FDR < 0.001) at one 
or more time points, 4515 DEGs were co-regulated by 
both stresses, accounting for 47.6% and 69.2% of the 
total DEGs, respectively. Of the co-regulated DEGs, 971 
DEGs had significant expression changes at all drought- 
and cold-treated time points, including 334 with 
up-regulation, 618 with down-regulation and 19 with 
complex regulation patterns. The results clearly showed 
an overlap between the DEGs in response to drought and 
cold. The top 100 co-regulated DEGs were listed in the 
Additional file 10. 

Of the 971 co-regulated DEGs, 69 DEGs (33 up- and 
36 down-regulated) were annotated as TPs, including 
the AP2/EREBP, NAG, WRKY, MYB, bHLH, C2H2, 
bZIP, GRAS, C3H, Dof, TCP, COL and GRF families. In 
addition, 46 co-regulated DEGs were annotated as pro- 
tein kinases and phosphatases, and 17 co-regulated 
DEGs as synthetases or components in plant hormone 
biosynthesis or signaling. Besides, many stress-related 
functional genes appeared in the co-up-regulated DEGs, 
such as those encoding the enzymes in flavonoid biosyn- 
thetic pathways, ripening related proteins, thaumatin-like 
or osmotin-like proteins, transporters, cell wall-associated 
proteins, and some types of PR proteins. 

The co-up-regulated DEGs were significantly aggregated 
into two GO biological processes: response to stress and 
secondary metabolite (mainly aromatic compounds or 
phenylpropanoids) biosynthesis or metabolism. How- 
ever, the co-down-regulated DEGs were significantly 
enriched in carboxylic acid metabolic process and fatty 
acid or lipid biosynthetic process. In the molecular 
function category, the co-up- or co-down-regulated 
DEGs were mainly related to catalytic activity or oxido- 
reductase activity. Unexpectedly, both the up- and 
down-regulated DEGs were not significantly enriched in 
any cellular component terms (Table 5). 

According to KEGG pathway analysis, the co-up- 
regulated DEGs were significantly enriched in eight 
pathways. Most of these pathways are completely or 
highly overlapped and actually belong to two pathways: 
microbial metabolism in diverse environments and bio- 
syntheses of secondary metabolites. In contrast, the co- 
down- regulated DEGs were mainly assigned to meta- 
bolic pathways and starch/sucrose metabolic pathway. 
The DEGs involved in carotenoid biosynthesis and 
bisphenol degradation were also significantly enriched 
(Table 6). Furthermore, about 30% of the co-regulated 



DEGs were unknowns or unannotated genes. They may 
play unique roles in response and adaptation to drought 
and cold environments of A. mongolicus. 

Differentially expressed genes specifically responding to 
drought or cold stress 

In total, 4972 and 2009 DEGs were specifically regulated in 
at least one time point during drought or cold stress, re- 
spectively. Of these DEGs, 1057 (445 up-regulated, 566 
down-regulated and 46 complicatedly regulated) or 1055 
(514 up-regulated, 541 down-regulated) had significant ex- 
pression changes across all three time points under drought 
or cold stress, respectively. The top 100 drought- or cold- 
specific DEGs were listed in the Additional files 11 and 12. 

Among the 1057 drought-specific DEGs, 91 DEGs en- 
code TPs, 66 encode protein kinases and phosphatases, 
and 24 are involved in hormone or cell signaling. Com- 
pared to the co-regulated DEGs, more TP gene fam- 
ilies were found in the drought-specific DEGs. In the 
1055 cold-specific DEGs, however, fewer genes encode 
regulatory proteins. Some distinct differences appeared 
between the drought- and cold-specific DEGs. Por in- 
stance, genes encoding trehalose-6-phosphate synthase, 
trehalose-6-phosphate phosphatase and early responsive 
to dehydration protein were only found in the DEGs 
specifically up-regulated in drought; while genes encod- 
ing acidic endochitinase, protochlorophyllide reductase, 
protein chloroplast import apparatus 2 and chloroplast 
omega-3 fatty acid desaturase appeared only in the 
DEGs specifically up-regulated in cold. In addition, 
about 23.7% of the drought-specific DEGs and 32.3% of 
the cold-specific DEGs were functional unknowns. 

The DEGs specifically up-regulated in drought were 
mostly enriched in the GO terms "Response to stress", 
"Response to stimulus", "Membrane" and "Antioxidant ac- 
tivity" (Additional file 13). However, the DEGs specifically 
up-regulated in cold were mostly enriched in the terms 
"Aminoglycan catabolic process", "Plastid", "Chloroplast" 
and "pattern binding" (Additional file 14). On the other 
hand, the DEGs specifically down-regulated in drought 
were more enriched in the GO terms "Small molecule bio- 
synthetic process", "Organic acid biosynthetic process" 
and "Iron ion binding" (Additional file 13); while the 
DEGs specifically down- regulated in cold were highly 
enriched in the GO terms "Lipid localization" and "Anti- 
oxidant activity" (Additional file 14). These results indi- 
cated that distinct molecular mechanisms were involved 
in A. mongolicus seedlings under drought compared to 
those in cold stress. 

Discussion 

A, mongolicus is an emerging model species for studying 
the response and tolerance of plants to extreme cold and 
drought environments. By using Illumina technology. 
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Table 5 Significantly enriched GO terms in the DEGs co-regulated by drought and cold^ 



GO term 



DEGs^ 



Unigenes^ 



Corrected P- value 



Co-up-regulated 
Biological process 

Response to stress 

Cellular aromatic compound metabolic process 

Defense response 

Secondary metabolic process 

Aromatic compound biosynthetic process 

Phenylpropanoid biosynthetic process 

Phenylpropanoid metabolic process 

Cellular component 

External encapsulating structure 

Chloroplast 

Cell wall 

Plant-type vacuole 
Chloroplast stroma 
Molecular function 
Catalytic activity 
Polysaccharide binding 
Pattern binding 
Carbon-nitrogen lyase activity 
Alkene binding 
Co-down-regulated 
Biological process 

Monocarboxylic acid metabolic process 

Fatty acid metabolic process 

Small molecule biosynthetic process 

Organic acid biosynthetic process 

Lipid biosynthetic process 

Fatty acid biosynthetic process 

Nucleobase transport 

Purine nucleoside transport 

Cellular component 

Intrinsic to membrane 

Cytoplasmic vesicle 

Vesicle 

Microbody part 
Peroxisomal part 
Molecular function 

Oxidoreductase activity 

Nucleobase transmembrane transporter activity 
Arsenate reductase activity 



36 (21.4%) 
15 (8.9%) 
15 (8.9%) 
1 1 (6.5%) 
10 (6.0%) 
9 (5.4%) 
9 (5.4%) 

9 (8.1%) 
7 (6.3%) 
6 (5.4%) 
3 (2.7%) 

2 (1.8%) 

135 (76.7%) 

3 (1.7%) 
3 (1.7%) 
3 (1.7%) 
2 (1.1%) 



27 (6.5%) 
24 (5.8%) 
23 (5.6%) 
22 (5.3%) 
19 (4.6%) 
17 (4.1%) 
9 (2.2%) 

5 (1.2%) 

27 (16.3%) 
17 (10.2%) 
17 (10.2%) 
2 (1.2%) 

2 (1 .2%) 

44 (18.6%) 

6 (2.5%) 

3 (1.3%) 



2565 (11.8%) 
418 (1.9%) 
580 (2.7%) 
232 (1.1%) 
154 (0.7%) 
147 (0.7%) 
200 (0.9%) 

857 (4.1%) 
614 (2.9%) 
329 (1.6%) 
109 (0.5%) 
56 (0.3%) 

15838 (64.3%) 
14 (0.1%) 
18 (0.1%) 
27 (0.1%) 
3 (0.0%) 



520 (2.4%) 
364 (1.7%) 
324 (1.5%) 
284 (1.3%) 
311 (1.4%) 
138 (0.6%) 
41 (0.2%) 
24 (0.1%) 

2357 (11.2%) 
1240 (5.9%) 
1256 (6.0%) 
33 (0.2%) 
33 (0.2%) 

2608 (10.6%) 
43 (0.2%) 
4 (0.0%) 



0.05017 
0.00018 
0.00935 
0.00040 
6.91 E-05 
0.00044 
0.00535 



1 

0.38977 
0.97678 
1 

0.03081 
0.01477 
0.03243 
0.11089 
0.01818 



0.00096 
5.22E-05 
2.53E-05 
1 .04E-05 
0.00316 
4.13E-07 
2.18E-05 
0.02775 

1 

0.99779 
1 
1 
1 

0.01967 
0.00043 
0.00044 



*The significantly enriched GO terms were determined using a corrected P-value<0.05. If no significantly enriched GO terms were identified, the top five terms 
were listed, ^number of DEGs assigned to certain GO term; "^number of all reference unigenes assigned to certain GO term. 
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Table 6 Significantly enriched KEGG pathways in the DEGs co-regulated by drought and cold^ 



Pathway 



DEGs^ 



Unigenes'^ 



Corrected P- value 



Pathway ID 



Co-up-regulated 

Microbial metabolism in diverse environments 29 (15.3%) 

Biosynthesis of secondary metabolites 28 (14.7%) 

Limonene and pinene degradation 17 (9.0%) 

Aminobenzoate degradation 17 (9.0%) 

Bisphenol degradation 9 (4.7%) 

Polycyclic aromatic hydrocarbon degradation 9 (4.7%) 

Chloroalkane and chloroalkene degradation 6 (3.2%) 

Phenylalanine metabolism 5 (2.6%) 
Co-down-regulated 

Metabolic pathways 55 (26.8%) 

Starch and sucrose metabolism 13 (6.3%) 

Carotenoid biosynthesis 7 (3.4%) 

Bisphenol degradation 7 (3.4%) 



2079 (8.8%) 
1995 (8.4%) 
1009 (4.3%) 
1030 (4.3%) 
221 (0.9%) 
248 (1.0%) 
187 (0.8%) 
108 (0.5%) 

4071 (17.2%) 
594 (2.5%) 
84 (0.4%) 
221 



0.03272 
0.03272 
0.03274 
0.03274 
0.00497 
0.00594 
0.03274 
0.03272 

0.00961 
0.04050 
0.00049 
0.04684 



ko01120 
koOlllO 
ko00903 
ko00627 
ko00363 
ko00624 
ko00625 
ko00360 

koOllOO 
ko00500 
ko00906 
ko00363 



*The significantly enriched KEGG pathways were determined using a corrected P-value < 0.05. ^number of DEGs assigned to certain KEGG pathway; number of all 
reference unlgenes assigned to certain KEGG pathway. 



9309 up-regulated and 23419 down-regulated DEGs were 
previously identified in the cold-acclimated A. mongolicus 
seedlings [12], Using an EST approach, 120 abiotic stress- 
responsive genes were identified in this species, and 82 
genes were confirmed to be cold- and/or drought- 
inducible by qRT-PCR [15]. In the present study, we per- 
formed comparative transcriptome profiling of A, mongo- 
licus in response to drought and cold and identified 
thousands of the DEGs regulated by each stress. In par- 
ticular, 971 co-regulated DEGs by both stresses and about 
1050 DEGs specifically regulated by a single stress were 
identified, suggesting common and distinct molecular 
mechanisms underlying the response to drought and cold 
environments. These results were similar to those found 
in Arabidopsis [7], rice [8] and Brachypodium [9] despite 
of some new information also observed in the present 
study. The co-regulated DEGs could represent a subset of 
basal drought and cold responsive genes in A, mongolicus 
and thus may play essential roles in its adaptation to ad- 
verse environments; while the specifically regulated DEGs 
might underline the distinct mechanisms of this species in 
response to different stresses. Functional enrichment 
analyses of these DEGs provided clues to the molecular 
bases of A, mongolicus combating drought and cold 
environments. 

Flavonoids may contribute greatly to the adaptation of ^. 
mongolicus to drought and cold environments 

Flavonoids constitute a large group of polyphenolic second- 
ary metabolites in plants, which have antioxidant activity 
and are of prime importance for plant defense against path- 
ogens and UV stress [16]. Recent studies have showed that 



flavonoids accumulated to high levels in high alpine and 
polar plants [17] and in A, mongolicus [18] during cold ac- 
climation, suggesting their important roles in cold adapta- 
tion of these plants. The biosyntheses of flavonoids require 
the enzymes involved in the general phenylpropanoid path- 
way and its flavonoid branch pathways, including phenyl- 
alanine ammonia-lyase (PAL), cinnamate-4-hydroxylase 
(C4H), 4-coumarate CoA ligase (4CL), chalcone synthase 
(CHS), chalcone isomerase (CHI), chalcone reductase 
(CHR), flavanone 3-hydroxylase (F3H), dihydroflavonol-4- 
reductase (DFR), leucoanthocyanidin dioxygenase (LDOX), 
anthocyanidin synthase (ANS), isoflavone synthase (IFS), 
2-hydroxyisoflavanone dehydratase (HID), isoflavone re- 
ductase (IFR), isoflavone-3'-hydroxylase (IF3H), flavonol 
synthase (FLS), O-methyltransferase (OMT) and UDP- 
glycosyltransferase (UGT), etc [19,20]. In the present study, 
genes encoding almost all of these enzymes were not only 
moderately or highly expressed in A, mongolicus seedlings 
in normal condition (average RPKM is 25.3), but also were 
significantly up-regulated by drought and cold stresses. GO 
terms and pathways relevant to flavonoids were highly 
enriched in the co-up-regulated DEGs identified in both 
drought and cold stresses (Tables 5 and 6). The phenylpro- 
panoid or flavonoid biosynthetic pathways are usually regu- 
lated by MYB and/or bHLH TFs [21,22]. Genes encoding 
these TFs were also included in the co-up-regulated DEGs 
identified in the present study. Although their exact func- 
tions need to be characterized, our results suggest that 
genes involved in flavonoid biosyntheses may be crucial in 
the adaptation of A, mongolicus to drought and cold envi- 
ronments. An increased production of reactive oxygen spe- 
cies is a common consequence in plants under abiotic 
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stresses, which usually damage cellular membranes and 
other cellular components resulting in oxidative stress and 
eventually cell death [23]. It is likely that a high level of fla- 
vonoids may protect A. mongolicus cells from oxidative 
stress originated from drought and cold, and hence they 
contribute to stress tolerance. 

Membrane proteins may play important roles in drought 
adaptation of A, mongolicus 

In plants, membrane transport and perception systems 
have essential roles in maintaining cellular homeostasis 
under adverse environmental stresses through cell-to-cell 
and/or organ-to-organ communication [24,25]. An in- 
creased expression of multiple transporter and channel 
protein genes in response to drought has been found in 
some plant species [7,24,26]. Over-expressing some of 
these genes, such as aquaporin OsPIP2-2 [27] and ABC- 
transporter AtABCG36 [28], enhanced plant tolerance to 
drought and/or salt stress. In the present study, the GO 
term "Membrane" was highly enriched in the DEGs up- 
regulated in drought (Additional files 5 and 13), and ap- 
proximately 70 membrane protein genes were found with 
up-regulation under drought stress. Most of these genes 
were annotated to encode transporters and channel pro- 
teins, such as those for ions, sugars, organic acids, water, 
nitrate and phosphate, suggesting their involvement in 
drought adaptation of A, mongolicus. 

Receptor-like protein kinases (RLKs), another group of 
membrane proteins, are important signaling components 
that mediate plant responses to developmental and envir- 
onmental stimuli [25]. Recent studies have shown that 
some RLKs, such as Arabidopsis RPKl, GHRl and 
AtCRK4S, responded to drought, ABA, salt and/or cold 
stresses and they were implicated in ABA and BR signaling 
pathways [29,30]. Genes encoding cysteine-rich receptor- 
like kinases and LRR-RLKs (leucine-rich receptor-like 
protein kinases) were highly induced in chilling-treated 
Chorispora hungeana [31]. In the present study, several 
RLK genes, such as those encoding LRR-RLKs, lectin-like 
receptor kinase, leucine-rich repeat receptor-like tyrosine- 
protein kinase and G-type lectin S -receptor-like serine/ 
threonine-protein kinase, were up-regulated in drought- 
stressed A. mongolicus seedlings. Although their functions 
and detailed molecular mechanisms are not clear, an in- 
creased transcription of RLK genes suggests they may play 
important roles in drought adaptation of A. mongolicus, 

Chloroplast plays a central role in cold adaptation of A 
mongolicus 

Physiological studies have confirmed that chloroplast 
plays a central role in plant cold acclimation and freezing 
tolerance, which acts as sensor of plants responding to 
cold and light stress as well as the target of many cold ac- 
climation processes [32]. In the present study, genes 



related to chloroplast were specifically enriched in the 
DEGs up-regulated in cold (Additional files 8, 9 and 14), 
covering almost all aspects of chloroplast. The largest 
group of these genes was involved in photosynthesis, en- 
coding chlorophyll a-b binding proteins, photosystem II 
and photosystem I components, cytochrome b6/f complex 
subunit VIII, thioredoxins, ribulose-1, 5-bisphosphate 
carboxylase/oxygenase large and small subunits, as well as 
protochlorophyllide reductase and Mg-chelatase. This re- 
sult reinforced previous findings in woody or shrub plants 
with strong cold tolerance, including A. mongolicus [15], 
rhododendron [33] and Populus euphratica [34], in which 
several photosynthesis-related genes were up-regulated in 
cold acclimation. However, it was different from those ob- 
served in herbaceous plants like Arabidopsis [35], Brachy- 
podium [9], barley [36] and tomato [37]. In these plants 
with moderate cold tolerance or cold sensitivity, signifi- 
cant suppression of photosynthesis-related genes was ob- 
served under cold stress. An up-regulated transcription of 
the photosynthesis-related genes might be an important 
mechanism underlying the maintenance of normal or 
higher photosynthetic efficiency of A, mongolicus during 
extreme cold winter. 

The second group of the up-regulated genes was anno- 
tated to encode early light-induced proteins (ELIPs), 
20 kDa chaperonin, chaperone protein DnaJ 11, protease 
Do-like 8 (DEGP8) and lipoxygenase 4 (LOX 4). ELIPs 
might protect chloroplasts from light-induced damage or 
photooxidative damage [38]. Chaperonins or chaperone 
proteins participate in proper folding of protein sub- 
strates, protein disaggregation and degradation [39]. 
DEGP8 located in the thylakoid lumen [40] is probably 
involved in the cleavage of photodamaged proteins. 
LOX 4 is involved in the production of jasmonic acid and 
functions in biotic stress response [41]. To our knowledge, 
these genes (except for chaperone protein DnaJ 11) were 
firstly reported with up-regulation in A.mongolicus under 
cold stress. 

Two other chloroplast-related genes up-regulated in 
cold stress encode chloroplast omega-3 desaturase and 
glycerol-3-phosphate acyltransferase. Increasing evidence 
has indicated that both genes participate in the desatur- 
ation of fatty acids in chloroplast lipids and play crucial 
roles in cold response and tolerance [42,43]. An in- 
creased expression of these genes might help to maintain 
suitable fluidity and stability of chloroplast membranes 
of A. mongolicus in cold environment. The remaining 
chloroplast-related genes up-regulated under cold stress 
may have roles in protein import, ion transport, ATP 
synthesis, starch degradation, protein synthesis and pro- 
cessing. All these genes could provide structural and 
physiological mechanisms for the protection of photo- 
synthetic apparatus, and thus for extreme tolerance of A. 
mongolicus in very cold winter. 
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Transcriptional regulatory network is involved in the 
adaptation of >A. mongoHcus to drought and cold 
environments 

TFs are crucial components in stress signal transduction 
pathways, which directly control the expression of specific 
sets of stress-responsive genes [44,45] . Of the A. mongoH- 
cus transcriptome, at least 160 TFs in 13 TP families were 
identified as DEGs during drought- and/or cold treat- 
ments, such as AP2/EREBP, NAG, WRKY and bHLH 
families. 

DREBs and ERFs (ethylene responsive factors) are two 
subfamilies in the AP2/EREBP superfamily [46] . It is well- 
known that DREBs play central roles in the regulation of 
cold- and drought-responsive genes and are important to 
the related stress tolerance in many plant species [47]. 
ERFs have been recently found to participate in the re- 
sponse of plants to abiotic stresses, and some ERFs can 
enhance drought, salt and/or freezing tolerance of their 
transgenic plants [48,49]. In the present study, 9 DERBs 
and 8 ERFs were up-regulated in the drought- and cold- 
stressed A, mongolicus seedlings respectively, highlighting 
the importance of these TFs in the adaptation of this spe- 
cies to drought and cold environments. 

The NAG family is another group of TFs involved in the 
tolerance of plants to abiotic stresses, especially to drought 
and high salinity [50,51]. Our study identified 8 putative 
NAG genes up-regulated under drought and cold. Similar 
results were obtained in the WRKY family, a well-known 
family involved in biotic-stress tolerance in plants. Recent 
evidence showed that the WRKY family also functions in 
response to abiotic stresses [52]. Our findings provide sup- 
ports to these results, suggesting an important role of the 
NAG and WRKY families in drought and cold tolerance. 

The MYB, bHLH and G2H2 families are also involved 
in plant abiotic stress tolerance [5,45,53]. In the present 
study, different members of these families were either 
up- or down- regulated under drought and cold treat- 
ments. This differed obviously from the three TF families 
discussed earlier which were almost up-regulated during 
the stress treatments. In summary, sophisticated tran- 
scriptional regulation could participate in the adaptation 
of A, mongolicus to adverse environments, and several 
families, especially the AP2/EREBP (DREBs and ERFs), 
NAG and WRKY families might be important in the 
processes. 

Conclusions 

A comprehensive transcriptome resource is generated 
for A, mongolicus and transcriptome profiles under 
drought and cold stresses are characterized. Gompara- 
tive transcriptome analysis identified many genes and 
pathways commonly or specifically regulated in the two 
stresses, implicating molecular basis of drought and cold 
tolerance in A, mongolicus, A number of drought- and 



cold-regulated genes were firstly identified in this study, 
including some unknown genes. These newly found 
genes could be important to A. mongolicus in combating 
extremely cold and drought environments in desert. Our 
data will facilitate further molecular studies on stress tol- 
erance of A. mongolicus and provide new insights into 
plant adaptation to harsh environments. 

Methods 

Plant stress treatment and sampling 

A. mongolicus seeds were collected from a single shrub 
grown in the desert area in Dengkou Gounty, Inner 
mongolia, Ghina. These seeds were surface-sterilized with 
Glorox solution and soaked in water for three days at 25°G 
and then sown in pots filled with sand. The pots were 
placed in a greenhouse at 25°G under long-day condition 
(16 h light/8 h dark cycle). The seedlings were watered 
every three days with 1/5 strength of Hoaglands solution. 
Three months after germination, the seedlings were 
treated as follows. 

For de novo transcriptome sequencing and assembly: 
the seedUngs were placed into a low temperature- 
programmable incubator (BD-PRX-IOOOA; -10°G to 
50°G) for cold treatment at 4°G to -8°G under a 16 h 
dim light/8 h dark cycle in a gradual cooling way: 4°G 
for 6 h, 2°G for 14 h, 0°G for 4 h, -2°G for 12 h, -5°G for 
8 h, and -8°G for 4 h, in total for 48 h; the cooling rate 
was 0.5°G per minute. This cold treatment followed by a 
3-day recovery period can lead to a sUght lethal damage 
to some leaves of the A, mongolicus seedUngs according 
to our preliminary experiments (unpublished data). The 
seedUngs were carefully removed from sand at 1 h (4°G), 
3 h (4°G), 8 h (4°G 6 h + 2°G 2 h), 16 h (4°G 6 h + 2°G 
10 h), 24 h (4°G 6 h + 2°G 14 h + 0°G 4 h) and 48 h (cold 
treatment finished) after the initiation of the cold treat- 
ment at 4°G respectively, and washed with pre-cooled 
water (in a refrigerator at 3°G) and collected. Drought 
treatment was performed as follows: the seedUngs were 
carefully removed from sand, washed with tap water, 
placed on filter paper and transferred in the same incu- 
bator for dehydration at 25°G under a 16 h dim light/8 h 
dark cycle. The seedlings were collected at the same 
time points as cold stress (after 1, 3, 8, 16, 24 and 48 h 
of dehydration). The untreated seedlings (0 h) were used 
as control (GK). 

For transcriptome profile analysis: the seedlings were 
placed into the same low temperature-programmable in- 
cubator for cold treatment in a gradual cooling way: 4°G 
for 4 h, 2°G for 12 h, 0°G for 2 h, -2°G for 2 h, -4°G for 
2 h, and -6°G for 2 h, in total for 24 h. The drought 
treatment was performed using dehydration on a filter 
paper as described earlier. The seedlings were sampled 
at 2, 8 and 24 h after the initiation of the cold or 
drought treatment respectively. The untreated seedUngs 
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(0 h) were used as control (CK). For each time point, 
four seedlings were sampled. To minimize the bias of 
random changes from a single sample, each time point 
was sampled from three independent stress treatments 
and the total RNA was pooled equally for further 
experiment. 

Young leaves were sampled from a single A. mongoli- 
cus plant grown in the desert area mentioned above in 
summer (late July) and winter (late December). All sam- 
ples from stress treatments, controls and the young 
leaves were immediately frozen in liquid nitrogen and 
then stored at -76°C for RNA extraction. 

RNA extraction and lllumina sequencing 

Total RNA was extracted using a modified Trizol method 
developed by our laboratory (Trizol agent plus 2% p-ME 
and 3 M KAc with pH 4.8). The quantity of RNAs was 
verified by an ultraviolet spectrometer {OD2eo/OD28o ra- 
tios with 1.84 to 2.04). RNAs were dissolved in DEPC- 
treated H2O. For de novo transcriptome sequencing and 
assembly, total RNAs extracted from each sample were 
pooled equally, and 20 (ig of the pooled RNA was used for 
the isolation of poly (A) + mRNA. For transcriptome pro- 
file analysis, poly (A) + mRNA was purified from 10 (ig 
total RNA from each cold- or drought- stressed seedlings 
and control, respectively. 

The mRNAs were isolated by using oligo (dT) magnetic 
beads and was interrupted into short fragments by adding 
fragmentation buffer. Taking these short fragments as 
templates, the first-strand cDNA was synthesized using 
random hexamer-primer and reverse transcriptase. The 
second-strand cDNA was synthesized using buffer, 
dNTPs, RNase H and DNA polymerase I. After the end 
reparation and ligation of adaptors, the products were 
amplified by PGR and purified using the QiaQuick PGR 
extraction kit (Qiagen, Valencia, GA) to create cDNA li- 
braries. The cDNA libraries were sequenced by Beijing 
Genomics Institute (BGI)-Shenzhen (Shenzhen, Ghina) on 
an lllumina HiSeq™ 2000 according to the manufacturers 
protocols. For de novo transcriptome assembly, the 90 PE 
(paired-end) strategy was adopted; while for transcriptome 
profile analysis, the 50 SE (single-end) strategy was used. 

Transcriptome assembly 

Raw reads from the library were filtered to remove the 
dirty raw reads if they contain adaptors, over 5% of un- 
known nucleotides or over 50% of low quality bases 
(quality value < 5). Glean reads were used for de novo as- 
sembly with Trinity, a short read assembling program re- 
leased recently [54]. Briefly, Trinity firstly combined 
reads with certain length of overlap to form contigs. 
These reads were then mapped back to contigs. Finally, 
Trinity connected the contigs and generated unigenes, 
Le, the sequences that cannot be extended on either end. 



Functional annotation and classification of the unigenes 

The unigenes were aligned to protein databases in the 
NGBI in the priority order of NR, Swiss-Prot, KEGG 
and GOG by BLASTX algorithm [55], with a cut-off E- 
value < le-5. Proteins with the highest sequence similar- 
ity were retrieved as functional annotations for a given 
unigene. For unigenes that did not align to any of the 
above databases, ESTScan software was used to predict 
their coding regions and corresponding amino acid se- 
quences [56]. To identify putative TFs, all unigenes were 
blasted against the plantTFDB (http://planttfdb.cbi.edu.cn) 
with a cut-off E-value < le-5, and the top hit for each query 
was applied. Unigenes with no significant hits in the 
plantTFDB but with annotations as TFs in the NR and 
Swiss-Prot databases (E-value < le-10) were manually 
selected and classified into the corresponding TF fam- 
ilies. To ascertain the sequence similarity of A, mongoli- 
cus to other plant species, total unigenes were searched 
by BLASTN against the PlantGDB-generated Unique Tran- 
scripts (http://www.plantgdb.org/prj/ESTGluster/) (similar- 
ity > 80% and E-value < le-5). 

The GO annotations of the unigenes were obtained by 
using Blast2GO program [57]. The WEGO software was 
applied to perform GO functional classification for all 
annotated unigenes and to show the distribution of gene 
functions of the species at macro level [58]. After align- 
ing the unigenes to the GOG database, the GOG func- 
tional annotations were acquired. The KEGG pathway 
annotations were obtained by aligning the unigenes to 
the KEGG database [59]. 

Transcriptome profile analyses 

The clean reads from each library were mapped to the 
assembled transcriptome sequences (used as reference 
genes) using SOAPaligner/SOAP2 [60]. Mismatches with 
no more than two bases were allowed in the alignment, 
and ambiguous reads with multi-position match were 
excluded. Gene expression levels were measured by the 
RPKM (the numbers of reads per kilobase of exon re- 
gion in a given gene per million mapped reads) values 
[61]. If there were two or more transcripts for a gene, 
the longest one was used to calculate its RPKM. To 
identif)^ genes responding to drought and cold, the num- 
ber of reads and RPKM for each coding region in the 
control and each stress -treated library were determined, 
and the ratios of RPKM value for each gene in six stress - 
treated libraries to that in control library were calcu- 
lated. Statistical significance of the differential expression 
value for each gene was determined according to Audic 
and Glaverie [62]. 

Genes were selected as DEGs according to the FDR 
values < 0.001 and at least a two-fold change (the absolute 
value of log2 Ratio > 1) in expression level between two li- 
braries. GO enrichment analysis was performed by 
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Additional file 3: Histogram of COG classification of the annotated 
unigenes. In total, 17100 annotated unigenes were assigned to 25 COG 
classes. 

Additional file 4: The top 100 DEGs commonly regulated by 
drought stress. The unigenes with absolute values of log2 Ratio > 1 and 
FDR < 0.001 were identified as DEGs. 

Additional file 5: The significantly enriched GO terms in the 
DEGs commonly regulated by drought stress. The significantly 
enriched GO terms were determined using a corrected P-value < 0.05. 
If no significantly enriched GO terms were identified, the top five terms 
were listed, ^number of DEGs assigned to certain GO term; ^number of 
all reference unigenes assigned to certain GO term. 

Additional file 6: The significantly enriched KEGG pathways in the 
DEGs commonly regulated by drought stress. The significantly 
enriched KEGG pathways were determined using a corrected P-value < 
0.05. ^number of DEGs assigned to certain KEGG pathway; ^number of all 
reference unigenes assigned to certain KEGG pathway. 

Additional file 7: The top 100 DEGs commonly regulated by cold 
stress. The unigenes with absolute values of log2 Ratio > 1 and FDR < 
0.001 were identified as DEGs. 

Additional file 8: The significantly enriched GO terms in the 
DEGs commonly regulated by cold stress. The significantly enriched 
GO terms were determined using a corrected P-value < 0.05. If no 
significantly enriched GO terms were identified, the top five terms were 
listed, ^number of DEGs assigned to certain GO term; ^number of all 
reference unigenes assigned to certain GO term. 

Additional file 9: The significantly enriched KEGG pathways in the 
DEGs commonly regulated by cold stress. The significantly enriched 
KEGG pathways were determined using a corrected P-value < 0.05. 
^number of DEGs assigned to certain KEGG pathway; ^number of all 
reference unigenes assigned to certain KEGG pathway. 

Additional file 10: The top 100 DEGs co-regulated by drought and 
cold. These DEGs had absolute values of log2 Ratio > 1 and FDR < 0.001 
at all drought- and cold-treated time points. 

Additional file 11: The top 100 DEGs specifically regulated by 
drought stress. These DEGs had absolute values of log2 Ratio > 1 and 
FDR < 0.001 at three drought-treated time points. 

Additional file 12: The top 100 DEGs specifically regulated by cold 
stress. These DEGs had absolute values of log2 Ratio > 1 and FDR < 0.001 
at three cold-treated time points. 

Additional file 13: The significantly enriched GO terms in the DEGs 
specifically regulated by drought stress. The significantly enriched 
GO terms were determined using a corrected P-value < 0.05. If no 
significantly enriched GO terms were identified, the top five terms were 
listed, ^number of DEGs assigned to certain GO term; ^number of all 
reference unigenes assigned to certain GO term. 

Additional file 14: The significantly enriched GO terms in the 
DEGs specifically regulated by cold stress. The significantly enriched 
GO terms were determined using a corrected P-value < 0.05. If no 
significantly enriched GO terms were identified, the top five terms were 
listed, ^number of DEGs assigned to certain GO term; ^number of all 
reference unigenes assigned to certain GO term. 

Additional file 15: Primer sequences used in the RT-qPCR assay. 

The primer pairs for 9 selected unigenes were shown. 



mapping all the DEGs to GO terms in the GO database 
(http://www.geneontology.org/) and the number of uni- 
genes included in each term was calculated. Hypergeo- 
metric test [63] was used to identify significantly enriched 
GO terms in DEGs compared with the assembled tran- 
scriptome of A. mongolicus using p-value < 0.05. Pathway 
enrichment was analyzed by mapping DEGs to the KEGG 
database (http://www.genome.jp/kegg). The number of 
unigenes involved in each pathway was calculated. By 
comparing with the assembled transcriptome of A, mogoli- 
cuSy the significantly enriched pathways were identified 
using p-value < 0.05. 

Validation of DEGs by RT-qPCR 

Eleven genes with different expression patterns were 
chosen for the validation. The template cDNAs were syn- 
thesized from 1.5 (ig of total RNAs (pre- treated with 
DNase I) from the same samples as those for the tran- 
scriptome profile analysis using M-MLV Reverse Tran- 
scriptase (Promega, USA). The 2 x SYBR-Green I RT-PCR 
Master Mix (Takara, Japan) was used as labeling agent, 
and a-tubulin of A, mongolicus was served as an internal 
reference gene. The reaction mixture (20 \\L) contained 
10 (iL 2x Master Mix, 10 (iM forward and reverse primers 
(Additional file 15) (0.4 (iL each) and 1 (iL template cDNA 
(diluted 10 folds with deionized water). The PGR program 
was as follows: 95°C for 30 sec; the first step: 40 cycles of 
95°C for 5 sec, 55 to 60° for 20 sec and 72°C for 20 sec; 
the second step: 1 cycle of 95°C for 15 sec, 55 to 60° for 
1 min, and 72°C for 15 sec. The reaction was performed 
on a Statagene Mx3000 PGR system. Three independent 
biological replicates were performed for each sample. The 
relative expression levels of the selected unigenes were 
calculated using the relative 2~^^^^ method [64]. Results 
represent mean standard deviations of the three experi- 
mental replicates. 

Availability of supporting data 

The data set supporting the results of this article is avail- 
able in the NCBI SRA (short read archive) repository 
under the accession number of SRS526744 (NCBI SRA, 
http://www.ncbi.nlm.nih.gov/sra?term= SRS526744). 

Additional files 



Abbreviations 

TF: Transcription factor; DEG: Differentially expressed gene; NR: Non- 
redundant; GO: Gene Ontology; COG: Clusters of Orthologous Groups; 
KEGG: Kyoto Encyclopedia of Genes and Genomes; ABA: Abscisic acid; 
BR: Brassinosteroid; FDR: False discovery rate; RT-qPCR: Real-time quantitative 
PGR; NAC: NAM, ATAF and CUC; DREB: Dehydration responsive element 
binding; HSP: Heat shock protein; GLP: Germin-like protein; TPP: Trehalose-6- 
phosphate phosphatase; CCD: Carotenoid cleavage dioxygenase; LEA: Late 
embryogenesis abundant; AP2/EREBP: APETALA2/ethylene-responsive 
element binding protein; MYB: Myeloblastosis oncogene; bHLH: Basic 
helix-loop-helix; C2H2: Cys2/His2-type zinc finger protein 2; bZIP: Basic- 



Additional file 1: Distribution of the annotated unigenes in 
each database, ^annotated or unannotated unigenes/total unigenes x 
100%; ^annotated unigenes in each database/total annotated unigenes x 
100%. 

Additional file 2: Species distribution of A. mongolicus sequences in 
the TIGR Plant Transcript Assemblies database. All unigenes were 
searched by BLASTN against the TIGR Plant Transcript Assemblies 
database. In total, 47486 unigenes had significant hits (similarity > 80% 
and E-value < le-5), and the sequence similarities to 12 species were 
shown in the figure. 
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domain leucine-zipper protein; C3H: Cys3/His-zinc finger protein; GRAS: GAT, 
RGA and SCR; TCP: TBI, CYC and PCF; COL: Constans-like; Dof: DNA binding 
with one finger; GRF: General regulatory factor; PR: Pathogenesis-related; 
PAL: Phenylalanine ammonia-lyase; C4H: Cinnamate-4-hydroxylase; 4CL: 
4-coumarate CoA ligase; CHS: Chalcone synthase; CHI: Chalcone isomerase; 
CHR: Chalcone reductase; F3H: Flavanone 3-hydroxylase; 
DFR: Dihydroflavonol-4-reductase; LDOX: Leucoanthocyanidin dioxygenase; 
ANS: Anthocyanidin synthase; IFS: Isoflavone synthase; HID: 2- 
Hydroxyisoflavanone dehydratase; IFR: Isoflavone reductase; IF3H: Isoflavone- 
3'-hydroxylase; FLS: Flavonol synthase; OMT: 0-methyltransferase; UGT: UDP- 
glycosyltransferase; RLK: Receptor-like protein kinase; LRR-RLK: Leucine-rich 
receptor-like protein kinase; BLIP: Early light-induced protein; DEGP8: Protease 
Do-like 8; LOX 4: Lipoxygenase; ERF: Ethylene responsive factor; RPKM: Reads 
per kilobase per million reads; SRA: Short read archive. 
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